% 这个octave脚本可用于简易演示二维扩散
% 冰淇淋融化模拟器（
% CC-BY 4

clc
clear

dx = 0.01;
dt = 0.00001;

D=2;
sigma = D*dt/(dx)^2;

[x y] = meshgrid(0:dx:2);
n = size(x, 1);

u0 = zeros(n,n);
u1 = zeros(n,n);
u2 = zeros(n,n);

for i = 1:n
  for j = 1:n
    if i > (0.5*n - 0.1*n) && i < (0.5*n + 0.1*n) && j > (0.5*n - 0.1*n) && j < (0.5*n + 0.1*n)
      u1(i,j)=1;
      u2(i,j)=1;
    end
  end
end

figure();
imgind=0;

diff_i = zeros(n-1);
diff_j = zeros(n-1);

for tick = 1:1000
  t(tick) = (tick-1)*dt;
  %  for i = 2:n-1
  %    for j = 2:n-1
  %      u2(i,j) = u1(i,j) + sigma*(u1(i-1,j)-2*u1(i,j) + u1(i+1,j)) + sigma*(u1(i,j-1)-2*u1(i,j) + u1(i,j+1));
  %    end
  %  end

  % 根据AI建议向量化
  diff_i = u1(1:n-2,2:n-1)-2*u1(2:n-1,2:n-1) + u1(3:n,2:n-1);
  diff_j = u1(2:n-1, 1:n-2)-2*u1(2:n-1, 2:n-1) + u1(2:n-1, 3:n);
  u2(2:n-1,2:n-1) = u1(2:n-1,2:n-1) + sigma*(diff_i+diff_j);

  u2(1,:) = zeros(1,n); %固定边界
  u2(n,:) = zeros(1,n);
  u2(:,1) = zeros(n,1);
  u2(:,n) = zeros(n,1);

	u0 = u1;
	u1 = u2;

  if mod(imgind,2) == 0
		clf
		hold on
		axis equal
    axis([0 2 0 2 -1 1])
    xlabel('x')
    ylabel('y')
    caxis([-1 1])

    view(0,90);
    z = u2;
    h = surf(x, y, z);
    set(h,'EdgeColor','none');
		drawnow
		pause(0.01)
  end
  imgind++;
end

